Overview

This document contains visualisation and analysis of data in Pak choi pollinators raw visits Brad correctedupdatedHeatherxlsx (002).xlsx which was previously cleaned.


This analysis was performed using R version 4.4.1 (2024-06-14 ucrt). The tidyverse suite of packages was used for data manipulation and visualisation. The emmeans package was used for model predictions.

Questions from Brad

  1. How does the presence and age of native SNH influence pollinator diversity (species richness and evenness and seed yield on exotic crop plants compared to bare fence lines on farms with the SNH (distance >400m from SNH) and bare-fence-lines on farms without SNH?

  2. Does farm type influence pollinator abundance and species composition and resulting seed yields?

  3. Does visitation differ between exotic or native species, or between bees and flies in relation to farm type or the presence of native plantings? Does this influence seed yield?

Data loaded

The cleaned data is loaded.


## tibble [5,865 × 11] (S3: tbl_df/tbl/data.frame)
##  $ farm_type                 : Factor w/ 2 levels "Arable","Dairy": 1 1 1 1 1 1 1 1 1 1 ...
##  $ farm_name                 : Factor w/ 19 levels "Ahuriri","Ashford_1",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ season                    : Factor w/ 3 levels "2021/22","2022/23",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ boundary_descriptor       : Factor w/ 6 levels "Control farm",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ boundary_age              : num [1:5865] NA NA NA NA NA NA NA NA NA NA ...
##  $ insect_key_number         : Factor w/ 51 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ new_insect_key_name       : Factor w/ 51 levels "Black bulb hover fly",..: 20 6 33 23 2 30 21 22 39 50 ...
##  $ sum_of_insects_100_flowers: num [1:5865] 2.58 0 0 0 0 ...
##  $ broad_insect_category     : Factor w/ 3 levels "Honey bees","Non-bees",..: 1 3 3 3 3 3 3 2 2 2 ...
##  $ origin                    : Factor w/ 3 levels "Native","Both",..: 3 3 3 1 3 1 1 3 2 3 ...
##  $ boundary_alt              : Factor w/ 4 levels "Control farm",..: 1 1 1 1 1 1 1 1 1 1 ...

Raw mean table

The data is summarised by taking the means for each insect, ignoring the treatments etc so we can get a feel for the numbers for each of the 51 insects.


new_insect_key_name broad_insect_category origin Mean_count_per_100_flowers
Lasioglossum bee Wild bees Native 7.378
Drone fly Non-bees Exotic 5.054
NZ orange hover fly Non-bees Native 4.866
Striped flesh fly Non-bees Exotic 1.586
NZ black hover fly Non-bees Native 0.987
Honey bee Honey bees Exotic 0.971
Seedcorn maggot fly Non-bees Exotic 0.870
Buff-tailed bumble bee Wild bees Exotic 0.621
Pollenia flies Non-bees Native 0.620
Other insects Non-bees Both 0.571
Ichneumonid wasps Non-bees Exotic 0.498
Bronze cluster flies Non-bees Exotic 0.469
Cabbage white butterfly Non-bees Exotic 0.436
Hylaeus bees Wild bees Native 0.241
Eleven spotted ladybird Non-bees Exotic 0.148
Other muscid flies Non-bees Both 0.113
European blue blow fly Non-bees Exotic 0.101
NZ blue hover fly Non-bees Native 0.098
Other parasitoid wasps Non-bees Both 0.097
Yellow admiral butterfly Non-bees Native 0.091
Other moths Non-bees Both 0.085
Globetail hover fly Non-bees Exotic 0.068
Grey triangle muscid fly Non-bees Native 0.064
Brown blow fly Non-bees Exotic 0.060
Streaktail hover fly Non-bees Native 0.059
European green blow fly Non-bees Exotic 0.055
Black Leioproctus bees Wild bees Exotic 0.051
Orange Leioproctus bee Wild bees Native 0.045
Magpie moth Non-bees Exotic 0.045
Small ginger blister fly Non-bees Exotic 0.042
Green soldier fly Non-bees Native 0.036
Black bulb hover fly Non-bees Exotic 0.034
Ginger blister fly Non-bees Native 0.031
Red admiral butterfly Non-bees Native 0.030
Three spotted fly Non-bees Exotic 0.025
Grey-black bristle flies Non-bees Native 0.025
Other hover flies Non-bees Both 0.023
Other wasps Non-bees Both 0.022
Common copper butterfly Non-bees Native 0.022
Other bristle flies Non-bees Both 0.016
Other bumblebees Wild bees Exotic 0.016
European tube wasp Non-bees Exotic 0.015
Blue muscid fly Non-bees Exotic 0.014
Other blow flies Non-bees Exotic 0.011
Flower longhorn beetle Non-bees Native 0.008
Scaptia horse flies Non-bees Native 0.006
NZ blue blow fly Non-bees Native 0.006
Other ladybirds Non-bees Exotic 0.005
Vespula wasps Non-bees Exotic 0.004
March fly Non-bees Native 0.004
Spotted hover fly Non-bees Native 0.004

Raw mean graphs

Broad insect category graph with ‘full’ boundary descriptors

The data is grouped by the broad insect category and the counts are summed, then the mean across the farm type, season and boundary descriptor is taken. These numbers are graphed.

Origin category graph with ‘full’ boundary descriptors

The data is grouped by the origin category and the counts are summed, then the mean across the farm type, season and boundary descriptor is taken. These numbers are graphed.

Broad insect category graph with ‘new’ boundary descriptors

The new boundary descriptor has the “1-3 year old” categories. collapsed into single category called New. Most of the time “One year old” is in 2021/22, “Two year old” is in 2022/23 and “Three year old” is in 2023/24.


In the first graph below, observations from the same farm are joined by a line.

Origin category graph with ‘new’ boundary descriptors

In the first graph below, observations from the same farm are joined by a line.


Remove data before modelling

We remove the six observations where 2023/24 is One/Two year old. This means that the category “new” just has


  • “One year old” in 2021/22
  • “Two year old” in 2022/23
  • “Three year old” in 2023/24


(If we don’t want to interpret anything to do with how “new” the plantings are, we can include this data. Let me know if this is the case.)

Modelling using the broad insect categories

Insects per 100 flowers


The raw means are given below


farm_type boundary_alt season broad_insect_category mean n
Arable Control farm 2021/22 Honey bees 1.49 2
Arable Control farm 2021/22 Non-bees 13.49 2
Arable Control farm 2021/22 Wild bees 2.02 2
Arable Control farm 2022/23 Honey bees 0.83 2
Arable Control farm 2022/23 Non-bees 2.52 2
Arable Control farm 2022/23 Wild bees 1.51 2
Arable Control farm 2023/24 Honey bees 0.31 5
Arable Control farm 2023/24 Non-bees 10.70 5
Arable Control farm 2023/24 Wild bees 2.88 5
Arable Bare fence 2021/22 Honey bees 0.35 3
Arable Bare fence 2021/22 Non-bees 21.34 3
Arable Bare fence 2021/22 Wild bees 0.94 3
Arable Bare fence 2022/23 Honey bees 0.33 3
Arable Bare fence 2022/23 Non-bees 9.69 3
Arable Bare fence 2022/23 Wild bees 3.96 3
Arable Bare fence 2023/24 Honey bees 0.92 4
Arable Bare fence 2023/24 Non-bees 7.45 4
Arable Bare fence 2023/24 Wild bees 6.33 4
Arable Old 2021/22 Honey bees 1.00 4
Arable Old 2021/22 Non-bees 27.76 4
Arable Old 2021/22 Wild bees 4.27 4
Arable Old 2022/23 Honey bees 1.86 4
Arable Old 2022/23 Non-bees 14.73 4
Arable Old 2022/23 Wild bees 5.73 4
Arable Old 2023/24 Honey bees 4.22 5
Arable Old 2023/24 Non-bees 19.77 5
Arable Old 2023/24 Wild bees 10.16 5
Arable New 2021/22 Honey bees 0.50 3
Arable New 2021/22 Non-bees 21.23 3
Arable New 2021/22 Wild bees 2.84 3
Arable New 2022/23 Honey bees 1.91 4
Arable New 2022/23 Non-bees 12.39 4
Arable New 2022/23 Wild bees 7.19 4
Arable New 2023/24 Honey bees 2.65 4
Arable New 2023/24 Non-bees 29.17 4
Arable New 2023/24 Wild bees 12.14 4
Dairy Control farm 2021/22 Honey bees 0.00 7
Dairy Control farm 2021/22 Non-bees 28.00 7
Dairy Control farm 2021/22 Wild bees 5.14 7
Dairy Control farm 2022/23 Honey bees 0.07 7
Dairy Control farm 2022/23 Non-bees 7.53 7
Dairy Control farm 2022/23 Wild bees 1.72 7
Dairy Control farm 2023/24 Honey bees 0.04 6
Dairy Control farm 2023/24 Non-bees 13.26 6
Dairy Control farm 2023/24 Wild bees 2.10 6
Dairy Bare fence 2021/22 Honey bees 0.38 7
Dairy Bare fence 2021/22 Non-bees 19.81 7
Dairy Bare fence 2021/22 Wild bees 6.35 7
Dairy Bare fence 2022/23 Honey bees 0.35 7
Dairy Bare fence 2022/23 Non-bees 4.92 7
Dairy Bare fence 2022/23 Wild bees 0.53 7
Dairy Bare fence 2023/24 Honey bees 2.06 3
Dairy Bare fence 2023/24 Non-bees 9.46 3
Dairy Bare fence 2023/24 Wild bees 2.05 3
Dairy Old 2021/22 Honey bees 0.71 8
Dairy Old 2021/22 Non-bees 25.17 8
Dairy Old 2021/22 Wild bees 20.70 8
Dairy Old 2022/23 Honey bees 0.92 8
Dairy Old 2022/23 Non-bees 8.48 8
Dairy Old 2022/23 Wild bees 17.21 8
Dairy Old 2023/24 Honey bees 0.94 8
Dairy Old 2023/24 Non-bees 24.26 8
Dairy Old 2023/24 Wild bees 22.34 8
Dairy New 2021/22 Honey bees 1.29 3
Dairy New 2021/22 Non-bees 46.42 3
Dairy New 2021/22 Wild bees 13.01 3
Dairy New 2022/23 Honey bees 0.07 3
Dairy New 2022/23 Non-bees 11.16 3
Dairy New 2022/23 Wild bees 4.29 3
Dairy New 2023/24 Honey bees 1.52 3
Dairy New 2023/24 Non-bees 28.82 3
Dairy New 2023/24 Wild bees 18.09 3


It is not possible to model the variability of the combination

Dairy, Control farm, 2021/22, Honey bees

because it has a mean of zero. This data is removed from the model.

We fit a zero-inflated negative binomial generalised linear mixed model. The log link function is used. The zero-inflated part of the model has just an intercept, the zero-inflation is the same for all the fixed-effect combinations. There is also a term in the model to account for over dispersion. The AIC was used to select that.

The GLMM part of the model has the four-way interaction of farm_type, broad_insect_category, boundary_alt and season There is a random effect for boundary_alt:farm_name.

The ANOVA table is


Chisq Df Pr(>Chisq)
farm_type 0.046 1 0.831
broad_insect_category 349.748 2 0.000
boundary_alt 47.376 3 0.000
season 81.309 2 0.000
farm_type:broad_insect_category 7.952 2 0.019
farm_type:boundary_alt 2.942 3 0.401
broad_insect_category:boundary_alt 34.891 6 0.000
farm_type:season 11.163 2 0.004
broad_insect_category:season 12.852 4 0.012
boundary_alt:season 9.856 6 0.131
farm_type:broad_insect_category:boundary_alt 23.814 6 0.001
farm_type:broad_insect_category:season 8.550 4 0.073
farm_type:boundary_alt:season 6.960 6 0.325
broad_insect_category:boundary_alt:season 11.935 12 0.451
farm_type:broad_insect_category:boundary_alt:season 11.004 11 0.443


Thankfully, the four-way interaction is not significant. We have one three-way interaction to consider

  • farm_type:broad_insect_category:boundary_alt

and two two-way interactions (with season in which is missing from the three-way)

  • farm_type:season
  • broad_insect_category:season


Number of taxa

The data is grouped by the broad insect category and the number of different insects with non-zero counts is calculated. These numbers are graphed. The honey bees and wild bees and put together.


We fit a negative binomial generalised linear mixed model. The log link function is used.

The GLMM part of the model has main effects for farm_type, broad_insect_category2, boundary_alt and season. There is a random effect for boundary_alt:farm_name.

The model did not converge with the four-way interaction, and the AIC was better with just main effects.

The ANOVA table is:


Chisq Df Pr(>Chisq)
farm_type 0.047 1 0.828
boundary_alt 52.917 3 0.000
season 22.583 2 0.000
broad_insect_category2 245.406 1 0.000


farm_type is not significant, but the other three factors are.


Modelling using the origin insect categories

Insects per 100 flowers

The raw means are given below


farm_type boundary_alt season origin mean n
Arable Control farm 2021/22 Native 9.70 2
Arable Control farm 2021/22 Both 0.99 2
Arable Control farm 2021/22 Exotic 6.32 2
Arable Control farm 2022/23 Native 1.40 2
Arable Control farm 2022/23 Both 0.00 2
Arable Control farm 2022/23 Exotic 3.46 2
Arable Control farm 2023/24 Native 6.15 5
Arable Control farm 2023/24 Both 0.24 5
Arable Control farm 2023/24 Exotic 7.49 5
Arable Bare fence 2021/22 Native 7.95 3
Arable Bare fence 2021/22 Both 0.00 3
Arable Bare fence 2021/22 Exotic 14.68 3
Arable Bare fence 2022/23 Native 8.11 3
Arable Bare fence 2022/23 Both 0.41 3
Arable Bare fence 2022/23 Exotic 5.45 3
Arable Bare fence 2023/24 Native 7.91 4
Arable Bare fence 2023/24 Both 0.19 4
Arable Bare fence 2023/24 Exotic 6.61 4
Arable Old 2021/22 Native 19.38 4
Arable Old 2021/22 Both 0.14 4
Arable Old 2021/22 Exotic 13.51 4
Arable Old 2022/23 Native 10.94 4
Arable Old 2022/23 Both 0.17 4
Arable Old 2022/23 Exotic 11.21 4
Arable Old 2023/24 Native 16.27 5
Arable Old 2023/24 Both 1.77 5
Arable Old 2023/24 Exotic 16.10 5
Arable New 2021/22 Native 10.30 3
Arable New 2021/22 Both 0.38 3
Arable New 2021/22 Exotic 13.89 3
Arable New 2022/23 Native 11.79 4
Arable New 2022/23 Both 0.49 4
Arable New 2022/23 Exotic 9.21 4
Arable New 2023/24 Native 16.89 4
Arable New 2023/24 Both 2.91 4
Arable New 2023/24 Exotic 24.17 4
Dairy Control farm 2021/22 Native 15.66 7
Dairy Control farm 2021/22 Both 0.28 7
Dairy Control farm 2021/22 Exotic 17.21 7
Dairy Control farm 2022/23 Native 3.92 7
Dairy Control farm 2022/23 Both 0.90 7
Dairy Control farm 2022/23 Exotic 4.49 7
Dairy Control farm 2023/24 Native 6.78 6
Dairy Control farm 2023/24 Both 1.49 6
Dairy Control farm 2023/24 Exotic 7.13 6
Dairy Bare fence 2021/22 Native 14.52 7
Dairy Bare fence 2021/22 Both 0.19 7
Dairy Bare fence 2021/22 Exotic 11.82 7
Dairy Bare fence 2022/23 Native 2.49 7
Dairy Bare fence 2022/23 Both 0.21 7
Dairy Bare fence 2022/23 Exotic 3.09 7
Dairy Bare fence 2023/24 Native 4.08 3
Dairy Bare fence 2023/24 Both 0.35 3
Dairy Bare fence 2023/24 Exotic 9.15 3
Dairy Old 2021/22 Native 28.92 8
Dairy Old 2021/22 Both 0.35 8
Dairy Old 2021/22 Exotic 17.31 8
Dairy Old 2022/23 Native 21.09 8
Dairy Old 2022/23 Both 0.29 8
Dairy Old 2022/23 Exotic 5.23 8
Dairy Old 2023/24 Native 27.65 8
Dairy Old 2023/24 Both 4.00 8
Dairy Old 2023/24 Exotic 15.89 8
Dairy New 2021/22 Native 50.56 3
Dairy New 2021/22 Both 0.53 3
Dairy New 2021/22 Exotic 9.64 3
Dairy New 2022/23 Native 7.26 3
Dairy New 2022/23 Both 2.32 3
Dairy New 2022/23 Exotic 5.94 3
Dairy New 2023/24 Native 23.70 3
Dairy New 2023/24 Both 1.99 3
Dairy New 2023/24 Exotic 22.74 3


It is not possible to model the variability of the combinations

  • Arable, Bare fence, Both, 2021/22
  • Arable, Control farm, Both, 2022/23


because the data has a mean of zero. We remove all the “Both” data from the model.


We fit a negative binomial generalised linear mixed model. The log link function is used. The model has the four-way interaction of farm_type, origin, boundary_alt and season. There is a random effect for boundary_alt:farm_name.

The ANOVA table is


Chisq Df Pr(>Chisq)
farm_type 0.424 1 0.515
origin 6.347 1 0.012
boundary_alt 49.261 3 0.000
season 75.967 2 0.000
farm_type:origin 4.566 1 0.033
farm_type:boundary_alt 3.441 3 0.328
origin:boundary_alt 14.221 3 0.003
farm_type:season 7.818 2 0.020
origin:season 2.149 2 0.341
boundary_alt:season 15.799 6 0.015
farm_type:origin:boundary_alt 4.772 3 0.189
farm_type:origin:season 1.455 2 0.483
farm_type:boundary_alt:season 4.511 6 0.608
origin:boundary_alt:season 4.090 6 0.665
farm_type:origin:boundary_alt:season 12.056 6 0.061


The four-way interaction is not significant. None of the three-way interactions are significant. We have four two-way interactions to consider

  • boundary_alt:season
  • farm_type:season
  • origin:boundary_alt
  • farm_type:origin


Number of taxa

The data is grouped by the origin insect category and the number of different insects with non-zero counts is calculated. We remove the “both” category because there too many zeros for modelling.

We fit a negative binomial generalised linear mixed model. The log link function is used.

The model has main effects for farm_type, broad_insect_category2, boundary_alt and season plus the two-way interaction between farm_type and season. There is a random effect for boundary_alt:farm_name.

The model has a parameterised dispersion structure using season.

The model did not converge with additional terms.

The ANOVA table is:


Chisq Df Pr(>Chisq)
farm_type 0.003 1 0.953
season 23.994 2 0.000
boundary_alt 50.457 3 0.000
origin 12.731 1 0.000
farm_type:season 8.656 2 0.013


The two-way interaction farm_type: season is significant, as are the main effects for origin and boundary_alt.


Modelling the numbers of five individual insects (those that were more abundant than honey bees)

The following were more abundant on average than honey bees: Drone fly, Lasioglossum bee, NZ black hover fly, NZ orange hover fly, Striped flesh fly.

In the five models below zero-inflation and dispersion were considered when fitting the models. The AIC and BIC were looked at in conjunction with simulateResiduals and other functions from the DHARMa package were used. In each case the model was “better” (or no worse) without additional terms.

Lasioglossum bees

Graphs of raw means and data used in model

Model

We fit a negative binomial generalised linear mixed model. The log link function is used. The model has the three-way interaction of farm_type, boundary_alt and season. There is a random effect for boundary_alt:farm_name.


The ANOVA table is


Chisq Df Pr(>Chisq)
farm_type 4.095 1 0.043
boundary_alt 51.212 3 0.000
season 4.612 2 0.100
farm_type:boundary_alt 9.105 3 0.028
farm_type:season 14.775 2 0.001
boundary_alt:season 4.746 6 0.577
farm_type:boundary_alt:season 9.001 6 0.174


There are two two-way interactions to consider

  • farm_type:season
  • farm_type:boundary_alt


Drone fly

Graphs of raw means and data used in model

Model

We fit a negative binomial generalised linear mixed model. The log link function is used. The model has the three-way interaction of farm_type, boundary_alt and season. There is a random effect for boundary_alt:farm_name.


The ANOVA table is


Chisq Df Pr(>Chisq)
farm_type 1.693 1 0.193
boundary_alt 4.003 3 0.261
season 29.422 2 0.000
farm_type:boundary_alt 5.964 3 0.113
farm_type:season 0.058 2 0.971
boundary_alt:season 8.180 6 0.225
farm_type:boundary_alt:season 10.923 6 0.091


The main effect for season is the only significant term.


NZ orange hover fly

Graphs of raw means and data used in model

Model

We fit a negative binomial generalised linear mixed model. The log link function is used. The model has the three-way interaction of farm_type, boundary_alt and season. There is a random effect for boundary_alt:farm_name.


The ANOVA table is


Chisq Df Pr(>Chisq)
farm_type 2.347 1 0.126
boundary_alt 1.195 3 0.754
season 36.430 2 0.000
farm_type:boundary_alt 0.915 3 0.822
farm_type:season 0.506 2 0.777
boundary_alt:season 11.071 6 0.086
farm_type:boundary_alt:season 4.278 6 0.639


The main effect for season is the only significant term.


Striped flesh fly

Graphs of raw means and data used in model

Model

We fit a negative binomial generalised linear mixed model. The log link function is used. The model has farm_type, boundary_alt and season plus all the two-interactions of these three terms. The model failed to converge with the three-way interaction so it was removed. There is a random effect for boundary_alt:farm_name.



The ANOVA table is


Chisq Df Pr(>Chisq)
farm_type 4.841 1 0.028
boundary_alt 4.549 3 0.208
season 34.651 2 0.000
farm_type:boundary_alt 0.844 3 0.839
farm_type:season 1.349 2 0.509
boundary_alt:season 9.202 6 0.163


The main effects for season and farm_type are significant.


NZ black hover fly

Graphs of raw means and data used in model

Model

We fit a negative binomial generalised linear mixed model. The log link function is used. The model has the main-effect for farm_type, boundary_alt and season. The interaction terms caused convergence problems (or were not significant). There is a random effect for boundary_alt:farm_name.


The ANOVA table is


Chisq Df Pr(>Chisq)
farm_type 10.649 1 0.001
boundary_alt 12.190 3 0.007
season 12.418 2 0.002


The main effects are all significant.